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Abstract 

Numerical techniques are described for discretization of velocity space in continuum 
kinetic calculations. An efficient spectral collocation method is developed for the speed 
coordinate - the radius in velocity space - employing a novel set of non-classical orthogo- 
nal polynomials. For problems in which Fokker-Planck collisions are included, a common 
situation in plasma physics, a procedure is detailed to accurately and efficiently treat 
the field term in the collision operator. When species with disparate masses are included 
C/3 ■ simultaneously, a careful extrapolation of the Rosenbluth potentials is performed. The 

(~| _ techniques are demonstrated in neoclassical calculations of the bootstrap current and 



plasma flows in a tokamak. 

Keywords: orthogonal polynomials, kinetic, Fokker-Planck, velocity space, phase 
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1. Introduction 

A ubiquitous situation in numerical kinetic calculations is that the distribution func- 
tion must be discretized in a manner allowing both accurate integration and accurate 
differentiation. Integration is needed because typically moments of the distribution func- 
tion such as density and velocity are of interest. Differentiation is needed due to both 
the collisionless terms in the kinetic equation and also for velocity-space diffusion in col- 
lisions. Spherical or cylindrical coordinates are natural for velocity space, meaning the 
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normalized spherical or cylindrical radius coordinate x lies in the semi-infinite domain 
[OjCXj). The distribution function often has a Maxwellian envelope, meaning it behaves 
as oc exp(— a;^) as x ^ oo. 

Many discretization schemes for the radial coordinate are possible and many have 
been explored in the literature [ll [3, [sj, each with advantages and disadvantages 
regarding the above requirements. A uniform grid allows modest accuracy at both in- 
tegration and differentiation using finite difference/ volume/element methods. To use a 
uniform grid, x may cither be mapped to a finite interval using a coordinate transforma- 
tion, or else the fact that the distribution function is exponentially small for x > 6 may 
be used to truncate the x grid above some Xmax- Alternatively, Gaussian abscissa permit 
extremely accurate integration but generally only low-order differentiation if finite differ- 
ence/volume/element methods are applied to the unequally spacedgrid. A Chebyshev 
grid permits both spectrally accurate differentiation and integration[y,Q|, but Chebyshev 
grids involve a high density of nodes near the cndpoints of a finite interval, and arc there- 
fore poorly suited for the semi-infinite domain of x and cx exp(— x^) dependence typical 
of distribution functions. Methods using Lagucrre or associated Laguerrc polynomials 
have also been used, but as we will show, these methods can perform surprisingly poorly 
due to a nonanalytic Jacobian in the coordinate transformation from energy to speed. 

In this work, we present a new approach for discrctizing velocity space. The approach 
permits both spectrally accurate integration and differentiation on the semi-infinite do- 
main [0, oo) for functions with Maxwellian-like cx exp(— x^) dependence for large normal- 
ized speed X. The method is collocation rather than modal in nature, i.e. the function is 
known on a set of grid points (abscissae) rather than being explicitly expanded in a set 
of modes. As such, the methods is well suited for nonlinear computations in addition to 
linear ones. 

Our method derives from a set of previously unexplored orthogonal polynomials. The 
semi-infinite integration domain in the orthogonality relation for these new polynomials 
is identical to that of the Lagucrre polynomials. However, the use of a different weight 
in place of e^^ results in polynomials with superior properties for the calculations 
of interest. 

In kinetic calculations for plasmas, it is often important to accurately include collisions 



in the kinetic equation, and the Fokker-Planck operator [sj is the best available description 
of collisions. This operator may be written in terms of "Rosenbhith potentials" , which 
are defined in terms of the distribution function through a pair of Poisson equations. 
A complication of the Fokker-Planck operator is that the Rosenbluth potentials vary as 
powers of x rather than as exp(— x'^) for large x, which means discretization schemes that 
work well for the distribution function may not work well for the potentials. Accurate 
numerical schemes for handling the Fokker-Planck operator in plasma computations have 
been a subject of great interest Q, ^, l^. Here we will develop an efficient approach 
to incorporating the Rosenbluth potentials in the kinetic equation, carefully accounting 
for their behavior at large x to maintain high precision even for coarse grid resolution. 

The new techniques we discuss are demonstrated in several applications. First, we 
compute the bootstrap current in a tokamak plasma. Secondly, we calculate the flows 
of multiple ion species in a tokamak. Both of these computations require the solution of 
drift-kinetic equations, equations in which both the collision operator and other kinetic 
terms appear. Using the new velocity discretization described here, we find that only 
4-6 grid points in x are required for the desired level of convergence. For comparison, 
previous numerical methods for the same problems [l^ Q, [l^ have typically used tens 
or hundreds of grid points or modes in x. Even for gyrokinetic simulations of plasma 
turbulence, which commonly use ^ 16 energy gridpoints, this new energy grid may result 
in a significant speedup. 

Despite the exceptional accuracy of the method reported here, finite difference and 
finite volume methods for differentiation do have the advantage that they can be de- 
signed to maintain conservation properties exactly [l^, which can be important in time- 
dependent simulations. A hybrid scheme can be used, with the new scheme used to 
determine the grid points and integration weights, with finite differencing used to evalu- 
ate derivatives maintaining exact conservation. 



2. Spectral collocation scheme for velocity space 

For a variety of problems in kinetic theory, cither with or without collisions, it is 

useful to represent the distribution function in either spherical or cylindrical coordinates 

in velocity space. The dimensional coordinates v (the spherical radius in velocity space) 
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or v±_ (the cylindrical radius in velocity space) then arise. Either coordinate may be 
normalized for numerical work by the thermal speed wth = y^2T/m where T is a typical 
temperature and m is the mass of the particle species. Wc define x ~ u/^th or a; = v±/vth 
as appropriate. 

For large x, distribution functions typically behave like a Maxwellian, cx exp(— x^). 
It is therefore natural to represent the distribution function as a sum of polynomials 
that are orthogonal according to the relevant weight and domain: 

P^(x)F„^(x)x^e-"' dx = M^SN,n (1) 

where k is any number greater than —1, and Af,^ is some normalization. Notice ((1} differs 

from the defining orthogonality relations for both the associated Lagucrre and Hermite 

polynomials, which arc, respectively, 

r°° r(v -\- m + \ ) 

L^{y)L^{y)y"^e-ydy=^—-^—^Sj,,„ (2) 

(i.e. polynomials in rather than x) and 

/oo 
HN{x)Hjx)e~''" dx^2''nlV^5N^n (3) 
-oo 

(different range of integration than ([Ij). Lagucrre polynomials are the special case of 
associated Lagucrre polynomials with m = 0. There are several reasons why it is prefer- 
able to use polynomials in speed x rather than polynomials in energy y ~ x'^ , i.c why the 
new polynomials are preferable to Lagucrre or associated Lagucrre polynomials. These 
reasons will be developed throughout the remainder of this section. As initial motivation, 
consider that the new polynomials can represent both even and odd powers of v ot v±, 
whereas the associated Lagucrre polynomials can represent only even powers. 

In our experience, the choice fc = tends to yield the fastest convergence for the 
problems wc consider in the following sections, so for the rest of this paper we consider 
the polynomials P„ = P^. It is straightforward to generalize all results and algorithms 
to the case of different k. 

The first few polynomials may be computed itcratively using the following Gram- 
Schmidt procedure (though this method turns out to be poorly conditioned when n is 
large). The polynomial Pn{x) has n+1 coefficients, which may be determined by imposing 

orthogonality with respect to Pq through P„_i and enforcing the normalization, for a 
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total of n + 1 constraints. The first few polynomials P„, normalized so the leading 
coefficient is 1, are thus 

P^{x) = 1 (4) 
1 



Pi(x) = 

P2{x) ^ X'- 



2 VTT 



7r-2 2(7r- 



^ , , n 37r-8 2 IO-Stt 16 - Stt 

P3(2:) = x^--^- + 777 " 



2V^(7r-3) 2(7r-3) 4V^(7r - 3) ' 
The same polynomials normalized so A/^ ~ 1 are 

M^) = ^ = 1.0623 (5) 

Pi(x) = 2.4921a;- 1.4060 

P2ix) = 4.2656a;2- 6.6228x+ 1.6037 

Psix) = 6.0027a;^- 17.0392x2 + 12.1931X- 1.7463. 

The appearances of y/n and tt in (|4]) highlight that the P„ polynomials do not seem to 
be simply related to any of the classical orthogonal polynomials. 

Figure [1] shows the first few normalized polynomials, multiplied by . Figure [2] 
shows the corresponding plot for the Laguerre polynomials. (Associated Laguerre poly- 
nomials look very similar.) These figures illustrate the first reason the new polynomials 
are preferable: they have structure localized more in the region where the particle den- 
sity is high {x < 2) compared to the L™(a;2) exp(— x^) modes. For example, the new 
polynomial modes have significant differences from each other in the range a; < 1 - so 
they can accurately resolve structure in this region where many particles reside - whereas 
the Laguerre modes are very similar to each other in the same range. Even with 10 La- 
guerre modes, no structure can be resolved within x < 0.5 since there are no nodes there, 
whereas the new modes permit ample resolution of this region. The same issue can be 
seen in figures [3] and IH which show the grids for the corresponding spectral collocation 
methods. In other words, the zeros of the Pn(x) arc plotted in Figure [31 and the square 
roots of the zeros of the L^n^\y) are plotted in Figure |4] to show the grid points in speed 
rather than energy. The Laguerre and associated Laguerre bases tends to "waste" more 

nodes in the region a; > 2 where the distribution function is nearly zero. 
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Figure 1: (Color online) The first 10 new polynomial modes P„(x) exp(— x^), normalized so AI^ = 1. 
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Figure 2: (Color online) The first 10 Laguerre modes L„(a;^) exp(— x^) (normalized). 

It should be mentioned that the abscissae can be scaled by a factor slightly less than 1 
without destroying the spectral convergence properties while achieving better resolution 
at small values of x. 

The clustering of nodes for the new polynomials near x = turns out to be a major 
advantage when multiple particle species are considered simultaneously, discussed in 
section O since high resolution is needed in this region when the species masses are 
disparate. 

As with any set of orthogonal polynomials, the set P„ is associated with a Gaussian 
integration scheme, with the abscissae corresponding to the zeros of the polynomials. 
These abscissae and Gaussian integration weights may be computed for any orthogo- 
nality relation using the "Stieltjes procedure" , an alternative method to construct the 
polynomials which is better conditioned numerically than the direct orthogonalization 
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Figure 3: (Color online) Zeros of the new polynomials P„{x), i.e. the speed grid for Gaussian integration 
and for the new spectral collocation method. 
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Figure 4: (Color online) Values of z = ^/y for zeros of the Laguorre polynomials L„{y) (dots) and 
associated Laguerre polynomials Ln^^\y) (crosses), i.e. the speed grids for Gaussian integration and 
for spectral collocation methods. 
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approach described above. The theory underlying the Stieltjes procedure is detailed 
in Refs. [l5, 16, 17|, with Fortran 77 subroutines implementing the algorithms in Ref. 
[isl l . Another practical implementation of the alg orithm may be found in Rcf. [2 : with 
documentation given in section 4.6 of Ref. 20|. In brief, the algorithm involves first 
iterativcly computing the coefficients aj and bj of the two-term recursion relation for the 
polynomials, 

Pj+iix) = (.T - aj)Pj{x) - bjPj^iix), (6) 
which may be done by evaluating 



!^Pj{x)e-^^dx 
S^P^{x)e-^'dx 
S^PlMe-^^dx 



(7) 



(8) 



The integrals in ©-([S]) are performed using adaptive quadrature, evaluating Pj{x) using 
the recursion relation. The resulting recursion coefficients are then used to form the 
tridiagonal Jacobi matrix: 



J = 



a.n-2 



\ 



\/K-i 



(9) 



As shown in the aforementioned references, the eigenvalues of this matrix give the ab- 
scissae. Furthermore, the first elements of the normalized eigenvectors may be squared 
and multiplied by dx — •\/7r/2 to obtain the integration weights. This algorithm 

is numerically well conditioned, and produces a grid with extremely accurate integration 
properties, as we will demonstrate shortly. 

The above procedure produces a grid with no point at x = 0. If a point is needed there 
to impose a boundary condition, i.e. if a Gauss-Radau grid is desired, the bottom-right 
element of the Jacobi matrix may be replaced by — 6„_iP„_2(0)/P„_i(0); as discussed 
in Section 6.2 of Ref. [17|. The abscissae and weights are then computed from the 



eigenvalues and eigenvectors as before. The resulting grid will still be extremely accurate 
for integration, but slightly less so than the grid with no point at a; = 0. 



As mentioned previously, for kinetic calculations we need not only an accurate inte- 
gration scheme for the grid, but also an accurate differentiation matrix. The spectral 
differentiation matrix for any set of nodes and any weig hting function may be computed 
using the poldif algorithm described in Refs. 2l|, |22|, |23|. The algorithm amounts to 



representing a function on the grid as a weighted sum of the polynomials, differentiating 
the polynomials, and evaluating the result on the grid. We apply the algorithm using 

2 

the weight , and using the abscissae calculated from Stieljes procedure above. 

To demonstrate the efficacy of the new collocation method, figures [5E] compare var- 
ious discretization methods for both integration and differentiation of several functions. 
We first choose the function a;^ cos(a; — 1/2) exp(— a;^) because 1) for large x it has a 
Maxwellian-likc envelope, 2) for small x it has the behavior typical of physical ve- 
locity moment integrands near a; = 0, and 3) due to the cos(a; — 1/2) factor it cannot 
be exactly represented by any of the polynomial approximation methods, and its Taylor 
series contains both even and odd powers of x. We measure the accuracy of differentia- 
tion by comparing the numerical derivative to the analytical derivative at whichever grid 
point is closest to 1. For the new method, results are shown for both the standard Gauss 
abscissae (no point at x = 0) and for Gauss- Radau points (including a; = 0). It can be 
seen that the new schemes perform substantially better in both integration and differ- 
entiation than other methods. The same results are reported in tabular form in Tables 
[T][21 We have performed a similar analysis to many functions besides that in figure Oa, 
and the behavior seen in[5]b-c is typical. 

The poor performance of both the Laguerre and associated Laguerre methods in this 
comparison may be understood as follows. These grids can allow very accurate integra- 
tion and differentiation, but only for functions that are analytic in y = y^, the argument 
of the polynomials in the orthogonality relation that defines them. Since -y/x is nonan- 
alytic at a; = 0, the accuracy of Laguerre and associated Laguerre methods fails for 
general functions of x. This issue has a practical implication for kinetic computations. 
In kinetic codes, it is common to want both the density moment (J cPv f cx dx x^ f = 
^ /q°° dy^f) and velocity moment (J d^vvf oc dxx^f = ^ dyyf) of the distri- 
bution function /. Even if / is analytic in y, the spectral accuracy of Gaussian integration 
is lost when computing the density moment using a Laguerre grid due to the ^^y factor. 



as noted in Refs. . This problem could be avoided using an associated Laguerre grid 
with m = 1/2, so the nonanalytic factor y/y is incorporated into the integration weight. 
However, then accuracy would be lost for the velocity integrand, since it includes an 
additional factor of By using polynomials in x rather than polynomials in y, good 
accuracy can be retained for both integrals simultaneously. 

The plots also include results for a uniform grid and Chebyshev grid covering the 
interval [0,5]. We choose reasonable cutoff since the model function is 

exponentially small (< 3x 10"^*^) for a; > 5. For the uniform grid, integration is performed 
using the trapezoid rule, and differentiation is performed using 5-point finite differencing. 

Also shown are results for several schemes in which x G [0, oo) is mapped to a new 
variable s in the finite interval [0, 1]. Wc considered the three maps s = — ln(l — x), 
X = tan(7rs/2), and x = s/(l — s), applying either a uniform or Chebyshev grid to s. 
Results are shown for the map which produced the most accurate results. The direct 
Chebyshev grids in x and s perform roughly equally well to each other, while the uniform 
grid in s tends to outperform the uniform grid in x. 

Results arc also plotted for the collocation scheme used in the gyrokinctic code GS2, 
detailed in Ref. Q]. The same grid is also used in the astrophysics code AstroGxIsl. In 
this scheme, the grid consists of ?i — m Gauss-Legendre points on the interval [0,2.5], 
together with m Gauss-Laguerre points on the interval [2.5, oo), where m ~ 1 for n < 12 
and m = 2 otherwise. The philosophy behind this approach is to only use the Gauss- 
Laguerre points away from the singularity at a; = 0. Indeed, the GS2 approach shows 
substantially better convergence of the integral than the Laguerrc-only or associated 
Laguerre methods. The actual GS2 fortran subroutines are used to generate both the 
abscissae and weights for the figures here. Convergence of the integral stagnates around 
lO"'' — 10~^ since no more than 2 points are ever used in the interval [2.5, oo); if m 
is increased beyond 2, the GS2 method can reach machine precision. The derivative is 
computed using a 3-point finite-difference stencil, just as is used for the second derivative 
in GS2. This low-order scheme was used in GS2 because it allowed exact numerical 
conservation of particles in velocity space, which higher-order methods generally do not 
allow. 

Figures [B][7] show interesting effects arising when the various integration and differen- 
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Figure 5: (Color online) A central result of this paper: for the function in (a), typical of the velocity 
dependence of a distribution function, the new polynomial spectral collocation method described here 
far outperforms other schemes for both (b) integration and (c) differentiation. In (b)-(c), n denotes the 
number of grid points. H 
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Figure 6: (Color online) Same as Figure[5l but for the function in (a). Missing points indicate zero error 
to machine precision. 
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Figure 7: (Color online) Same as Figure[5l but for the function in (a). Missing points indicate zero error 
to machine precision. 
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Figure 8: (Color online) Same as Figure [S] but for the function in (a). 
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Figure 9: (Color online) Same as Figure [S] but for the function in (a). 
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Discretization scheme n = A n = 8 n = 16 
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Table 1; Error in the integral of cos(x — l/2)e ^ for various numbers of grid points n. 
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Table 2: Error in the derivative of eos(x — l/2)e ^ for various numbers of grid points n, reported at 
the grid point closest to a; = 1. 
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tiation schemes are applied to x°'e~^ for various integers a. The new polynomial scheme 
gives answers correct to machine precision both for integration (when a < 2n) and for 
differentiation (when a < n). The Lagucrre or associated Laguerrc grids with integer m 
can get the derivative correct to machine precision when a is even, or they can compute 
the integral correct to machine precision when a is odd, but they perform poorly in the 
opposite cases. This behavior follows because the Lagucrre method is designed to ex- 
actly evaluate dyy^ e^^ = 2 dx x^^^^e^^ and to exactly differentiate y-'e"^ for 
integers but they will perform poorly when the nonanalytic factor is included in 
the integral or derivative. Conversely, associated Lagucrre grids with half-integer m can 
compute the integral correct to machine precision when a is even, or they can compute 
the derivative correct to machine precision when a is odd. but they perform rather poorly 
in the opposite cases. The new polynomials get all cases correct to machine precision 
for both even and odd a. In figure ([7]).b, the GS2 grid can reach machine precision for 
integration of this special integrand because on the interval [2.5, oo), the function can be 
exactly integrated by the Laguerre grid, even with just 2 points. Convergence is therefore 
limited only by convergence of the Legendrc grid on the interval [0, 2.5]. 

The integration and differentiation properties of the Lagucrre and associated Laguerre 

2 

grids on x^e"^ for even and odd a extend in part to all even and odd functions of x. 
The latter is illustrated in figure [51 which shows the performance of the various grids for 
the odd function j{x) = xJo{x) exp(— x^), where Jq is a Bessel function. (We choose j{x) 
since functions like it arise in gyrokinetics.) As xJo{x) may be expanded in a Taylor series 
about X ~ with only odd powers of x, j{x) is integrated accurately on a Laguerre grid, 
but not differentiated accurately on this grid. Similarly, j{x) is differentiated accurately 
on an associated Laguerre grid with half-integer m, but it is not integrated accurately 
on this grid. If the x preceding Jq is replaced by an even power of x, the behaviors of 
integer and half-integer associated Laguerre grids are reversed. However, as shown in 
figure |S1 the new polynomial grids allow very accurate integration and differentiation 
simultaneously. 

We do find that for some functions with small-scale oscillatory structure, the new 
polynomials may not be the most accurate discretization option. Such a function is 
illustrated in figure [SI For this function, all grid types show slower convergence of both 
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the integral and derivative compared to previous functions. (Note the expanded n axis 
in figure The Chebyshev grid in the remapped variable performs best, and the 

uniform remapped grid also performs well. The superior performance of these grids can 
be understood in part from their very high density of points near a; = 0, higher even than 
the new polynomials, helping to resolve the fine structure in this region. For problems 
in which fine structure is expected in the x coordinate, such as gyrokinetic turbulence, it 
may be useful to build codes that can easily switch between several types of grids, since 
it may not be clear a priori which will actually yield faster convergence in practice for 
the specific problem. 



3. Rosenbluth potentials 

Wc now consider issues related to the linearized Fokker-Planck collision operator [sj, 
which arises in many kinetic theory problems including neoclassical and turbulence cal- 
culations. For many problems of interest, the distribution function fa of each species a is 
nearly MaxwcUian: fa = /ao + fai with fai < /ao where faa Uaiir^/'^vl)-'^ exp(-a;^), 
ria is the density, Va ~ \/2Ta/ma is the thermal speed, and Xa = v/va- In this case, 
the operator Cab in the kinetic equation of species a due to collisions with species h may 
be linearized about the Maxwellians, and it is given by the sum of the test parti cle part 
Cabifaii fbo} and the field particle part Cab{faO, fbi}- The test particle part isji" 
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Cabifal, fbo} — ^ablT-b 

where Tab 



eTi{xb) 



1 d 



dv 



C{fal} + 'i^^xlfal + 



rUb 



dfai 
dv 



(10) 



crf(2;h) — 2n ^^^XbC 



AnZjZ'^e-^ hi A/ml, * 
drasekhar function, erf (a:;,) = 2n^^/^ J^'' e^"" dw is the error function. 



/{2xl) is the Chan- 



(11) 



is the Lorentz operator, ^ = v^^/v, v^^ = v ■ B/\B\ is the component of the velocity along 
the magnetic field, and the gyrophasc (p is the cylindrical angle in velocity space. The 
field particle part is 

"2i;2 d'^G 2v (ma 



Cab{faO, fbl} — ^abfaO 



V* dv"^ 



- 1 



rrib 



dv vl 



47r fb\ 

rub 



(12) 
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where H and G are the perturbations to the Rosenbluth potentials of species 6, defined 
by 



(13) 
(14) 



with the velocity-space Laplacian = \{d / dv)v'^ (d / dv) + 2£] . (The concise form 
of the field operator in ([T^ is derived in Eq. (7) of Ref. 125[-) For this paper we fo- 
cus on drift-kinetic calculations, so wc need not retain the gyrokinetic modifications to 
the test particle 



26[ and field particle 



25| operators that arise when the collision oper- 



ator is gyroaveragcd at fixed guiding center position. We also specialize to the usual 
case in which only the gyrophase-indcpendcnt distribution function is considered. (The 
gyrophase-dependent part of the distribution function is known analytically.) 

For the moment, consider a species affected by collisions only with itself (such as the 
ions in a pure plasma), so we may write the linearized collision operator as Cij/ai} — 
Caaifai, fan} + Caa{fan,fai}-, and WC drop the spccies subscripts to simplify notation. 
Suppose also we are solving a linear kinetic equation of the form 



(15) 



where D is some operator describing the kinetic terms and R is some inhomogcneous 
term. As we do not know the perturbed Rosenbluth potentials associated with /i, we 
may think of (|15p as having the block matrix structure 

/ 



Mil Mi2 Mi3 
A/21 A/22 
M32 M33 



\ 






1 R 


\ 




H 









J 


) 




I 





(16) 



The first row in this system corresponds to the original kinetic equation (jlSp . The second 



row corresponds to the Poisson equation (jl3|) , and the third row similarly corresponds to 

([H)) . The block structure was pointed out previously in Ref. There is no A/23 

or A/31 element in (|T5)) since G does not appear in (|T5)) and fi does not appear in 

We may identify Afn with the operator D in (fTSj) . plus the test-particle collision terms, 

and also the last term in (|12p (which does not involve H or G.) Also, A/12 consists of 

the H and dH/dv terms in A/13 is the G term in (dH), A/21 = 47r, A/22 = with 
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appropriate boundary conditions for H, AI^2 = —2, and M33 = with appropriate 
boundary conditions for G. 

For several reasons, it is beneficial to discretize the ^ dependence of H and G us- 
ing Legendre polynomials pi,{^): H{z,^) = Y.e He{z)pi{^) and G(z,0 = Y.eGi{z)pi{^) 
where z represents any other independent variables. The reasons stem from the fact that 
the Legendre polynomials arc cigcnfunctions of V|, : 



V ^1 

dv dv v'^ 



(17) 



for any W{v). First, the Legendre representation is extremely efficient, since V|, ~ i'^ 
so ^ 1/^^, and therefore very few modes need to be retained. Second, the Leg- 
endre discretization allows efficient treatment of the boundary at large v, which can 
be understood as follows. The distribution function will be as small as machine pre- 
cision for a; > 6, so it is inefficient to locate any grid points beyond this value of x. 

2 

However, H and G vary not as but rather as powers of x for large x (as we will 
show in a moment), so they arc not small on the scale of machine precision even for 
X > 6. In fact, \G\ generally increases to infinity with x, so it is invalid to impose any 
Dirichlet or Neumann boundary condition on G at large finite x. However, using the 
Legendre expansion above, for x > .XMax = 4—6, the defining equation for H becomes 
V%H = 0, and so Hi ~ AiV^^^^^^ + Bev^. The Bi ^ solutions are not physical, and 
so the Robin boundary condition vdHi/dv + (£ + l)H( ~ may be applied at XMax 
to ensure the correct behavior cx v~^^'^^\ For the G potential, the defining equa- 
tion V^G = 2H admits four linearly independent solutions. Two are the homogeneous 
solutions to V^G = 0, and the other two arc particular solutions, varying as the homo- 
geneous solutions times v^. Thus, G^ = G^w"'^"'"^-' + D(V^ + Eev^~^ + F^v^^"^ . As can be 
seen in Eq. (45) of Ref. [s], the 7^ and Fi ^ solutions are unphysical. To retain 
both of the remaining two solutions, the boundary condition must be second order, and 
can be found as follows. Writing d^Gi/dv^ + (iv dGi/dv + ct G^ = 0, the requirements 
that Gi (X and Gi oc v^~^ be solutions yield two equations for (/3, a). Thus, the 

condition v'^ d^Gi/dv^ + {2£+ l)vdGe/dv+ (£^ — 1)G£ = is the one we impose at XMax- 
Legendre polynomials tend to be efficient for representing fi as well, although for 
some applications such as bounce-averaged codes, other representations of the pitch- 



angle dependence of fi may be necessary and/or preferable 
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13|, l28|. For simplicity. 



we assume for the rest of this discussion that the ^ dependence of /i is discretized in the 
same Legendre representation as the potentials. If a different discretization is used for 
the 5 dependence of /i, A/12, and M21 must contain the appropriate matrices for 

mapping between the representations. 

It is much faster to invert the M22 and M33 blocks of (|T6)) than to invert Mu, for 
two reasons. First, it typically takes many more Lcgcndrc modes to represent /i than 
to represent the potentials. Secondly, the kinetic operator D included in Mu generally 
involves coupling between real space and velocity space coordinates, whereas M22 and 
M33 are independent of real space. In fact, since M22 and Af33 are diagonal in the 
Legendre index, the only coordinate with off-diagonal terms in these blocks is x. Thus, 
it can be convenient to formally solve the system (|16p before discrctizing, yielding the 
form 

(Mil - [M12 - MuM^^Hd32] M^^'M2i) h = R- (18) 

The fact that M22 and 1/33 are very sparse and fast to invert means that extremely 
high resolution in x can be used for H and G for almost no extra cost. This result 
is fortuitous, because the spectral collocation method described in the previous section 
cannot be used for H or G. This is because, as discussed earlier, the potentials behave 
as powers of x rather than as cxp(— x^) for large x. Consequently, H and G may be 
discretized using a very high resolution uniform grid in x with finite-difference differen- 
tiation derivatives used in Af22, ^33, and Afi3. Or, a high resolution Chcbyshcv grid in 
X could be used for the potentials, with spectral differentiation matrices used in Af22, 
M33, and Afi3. (For applications in following sections, we use a uniform grid.) In either 
case, the high resolution grid no longer appears once the matrix multiplications in (jl8p 
are performed, so the high resolution docs not affect the time to solve the system 

A complication introduced by the use of different x grids for /i and the potentials 
is that functions must be mapped from one grid to the other in A/12, A/13; and A/21. 
Mapping from the polynomial grid to the uniform grid can be done with spectral accuracy 
by evaluating the new polynomials on the uniform grid . An efficient algorithm for this 



purpose is the polint subroutine from Rcfs. 
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22I [23[ . To map from the uniform grid 



to the polynomial grid, we may use cubic interpolation. The limited accuracy of this 
interpolation method is not a problem because the grid resolution for the potentials can 
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be extremely high at httle cost. Or if the potentials are discretized on a Chebyshev grid, 
the polint subroutine may again be used to interpolate with spectral accuracy. To form 
Mi2, Mi3, and M21, both the "/i ^potentials" and "potentials—)- /i" interpolations 
must be expressed as explicit matrices. As a specific example of these concepts, consider 
that the matrix A/13, which corresponds to the G term in (|12l) . may be constructed as 
A/13 = ErA where E — Tabfaoi'^v'^ /v^) is a diagonal operator on the /i grid, r is the 
matrix representing cubic interpolation from the potential grid to /i grid, and A is the 
finite-difference matrix representation of d^/dx^ on the potential grid. 



4. Application 1: Bootstrap current 

For a concrete application of all the aforementioned numerical techniques, we consider 
a calculation of the bootstrap current, a current j driven along the magnetic field B due 
to radial gradients of density and temperature in a toroidal plasma. In an axisymmetric 
(tokamak) plasma consisting of electrons and a single ion species, the bootstrap current 
29,13,13 



has the form 



{j ■ B) = -dp. 



I / dpc dpi\ £32 dTc , £34" dTi 
^31 — -TT + -T7 + — rr 



(19) 



Pc \ dip dip J To dip ZTc dip 
using notation consistent with Ref. [sJl. In (jTH)) . c is the speed of light, / = RB(^ is 
the major radius R times the toroidal magnetic field Bq, Z is the ion charge, Tj and 
are the ion and electron temperatures, pi and p^ are the ion and electron pressures, and 
2'Kip is the poloidal magnetic flux. The quantity a is a number obtained by solving the 
ion kinetic equation, and the quantities £31, £32, and £34 are moments of the electron 
distribution function for various driving terms, with each distribution obtained by solving 
the electron kinetic equation. Specifically, 



£31 - -47r-i/2 (b d^£ dx x^fl^' 



(20) 



where f^^ is the solution of the normahzcd electron "drift-kinetic" equation 



^ dQ 2B \de J qR^b-Ve 2 de ^ ' 

where x = Xc, and the independent variables are x, ^, and the poloidal angle 0. Here, 

b = B/B, B = \B\, B — B/Bq, Bq — I/Ro, Ro is a typical major radius, q is the safety 
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factor, is the poloidal magnetic flux, i>' = VcqRol'i'c is a dimensionless measure of 
coUisionality, Vc ~ 4\/27rnce^ In A/(3-y/m^To^^^) is the electron coUision frequency, Vc = 



■\/2To/?Tic is the electron thermal speed, and = 1^^^ [Ccc + Cd] ■ Angle brackets denote 
a flux surface average: 

{Y) = {V'y^ / d0Y/B ■ Ve (22) 
Jo 

for any quantity Y where 

= ^ d9/B -^6 = j die/ Be, (23) 

d(.0 is the poloidal length element, V' — dV/dip, and 2TTV{tp) is the volume enclosed by 
the flux surface. 



By expanding Cd in ^Jm^Jra-i ~ 1/60 for deuterium, it can be shown that the electron- 



ion collision operator is dominated 
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2J] by the £ term in the test particle operator 



([To)) . Then since erf(x) — ^{x) — > 1 for .t ^ 1, the collision operator to use in ((2T|) is 
C'{/if^} = i^-Mac{./e^\/co} + Cec{/co,/if^}+reiniz;-3£{/W|]^ As the non-Maxwcllian 

ion distribution function /ii and perturbed ion Rosenbluth potentials are then no longer 
present in ((2T|) . we may solve the single-species electron equation ([2T|) without coupling 
to the ion equation. 

The other coefficients in the bootstrap current are calculated in a similar manner: 



-32 



4n-'/'(B di dxx^U^J^^) (24) 



where f'^"^ is the solution of (j21l) with an extra factor — 5/2 multiplying the right-hand 
side, and 

C,,^A^-^'^{B^y' (b j'^d^J^^dxxHf^^'^'j (25) 

where f^'['^ is the solution of (pi]) with the right-hand side replaced by (3^^ — l)(a;^/2)e~^ dB/dO. 

We now demonstrate the solution of ((2T|) and the corresponding equations for £32 and 
£34 using the velocity discretizations described in the preceding sections. Some model 
must be supplied for the 9 dependence of the dimensionless quantities B and qRob ■ \/6. 
For this example we use the Miller equilibrium model [sfj detailed in the appendix. 

The coordinate may be discretized in several ways. In decreasing level of spar- 
sity, some options are finite differencing, a Fourier modal expansion, or Fourier spectral 
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collocation[7|. (The modal approach need not be fully dense in 9 since many coupling 
terms will be exponentially small and can therefore be set to zero.) For results shown 
here we use the modal approach. The finite-difference and collocation approaches lead 
to simpler code, since multiplications of the unknown fields by functions of 9 may then 
be done pointwise, whereas the modal approach requires matrix multiplication. 

If a Gaussian grid is used there is no point at x = 0, or if the Gauss-Radau grid is 
used, there is a grid point at a; = 0. In the former case there is no boundary condition 
imposed at x — 0, and in the latter case, the boundary condition df^^ /dx = at a; = 
is imposed for the t ~ Q Lcgendre mode. We find results arc nearly identical for the two 
approaches. Results shown below use the grid with no point at a; = 0. No boundary 
condition needs to be applied f^Ci at x = oo, since the behavior is automatically 
enforced by the choice of basis functions. 

The reduced system (|18p may be solved efficiently using a sparse direct solver, or when 
the resolution is sufficiently high, faster solution is typically possible using an iterative 
Krylov-space solver. For the iterative solvers, preconditioning is essential. We find an 
effective prcconditioner can be obtained by eliminating all off-diagonal blocks in the x 
coordinate. Calculations shown here are performed using the BICGstab(l) algorithm 

H- 

Figure fTOl shows the convergence of the code outputs £31, £32, and £34 with respect 
to the various numerical parameters, using the new polynomial grid with no point at 
X = 0. For this run we choose the coUisionality = v' A^/'^ to be 0.01, where A = 3.17 
is the aspect ratio. The excellent convergence with respect to all coordinates is evident. 
For the scan of XMax, Ny is varied proportional to .TMax to maintain constant resolution. 
The very rapid convergence with the grid size in x is enabled by the efficient spectral 
methods described in the preceding sections. Convergence to two digits of precision 
is plenty for comparison with experiment or for experimental design, and the code is 
converged beyond this level with only 4-5 grid points in x. The code execution time is 
nearly independent of Ny, the resolution of the uniform x grid for the potentials, for the 
reasons discussed in section[21 Fortran /Petsc jssl Is^ and Matlab versions of the code run 
in comparable time, since the rate-limiting step is iCZ-factorization of the prcconditioner, 
which is highly optimized in Matlab. Execution times reported in Figure [TU] are for the 
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Figure 10: (Color online) Convergence of the bootstrap current coefRcicnts £31, £32, and C34 in II19II 
from the single-species drift-kinetic code with respect to numerical parameters: Nx = number of grid 
points in x for the polynomial grid, Ny = number of grid points on the uniform x grid for the Rosenbluth 
potentials, A'^^ = number of Legendre modes used to represent the distribution function, NL = number of 
Legendre modes used to represent the potentials, N6 = number of Fourier modes in 8, and xjviax = cutoff 
of grid for the potentials. Except for quantities being scanned, the numerical parameters are Nx = 6, 
Ny = 20, Ni; = 100, NL = 4, N9 = 25, and XMax = 5. The convergence at a very small Nx enabled 
by the new discretization scheme described here is a central result of this paper. The coUisionality is 
= 0.01. 
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Matlab version running on a Dell Precision laptop with Intel Core 17-2860 2.50 GHz CPU 
and 16 GB memory. 

Notice how few Legendre modes 2) need to be retained in the potentials for 
good convergence, whereas many more Legendre modes (~ 100) are required for the 
distribution function. The fact that much higher resolution is needed in the ^ and 9 
coordinates than in the x coordinate for i^* ^ 1 is a result of a boundary layer known 
to exist in the solution[35|, illustrated in Figure 5 of Ref. 36|. This sharp features lies 
along the curve in the (0,^) plane corresponding to the boundary between trapped and 
circulating particle orbits. The boundary layer widens as coUisionality is increased, so 
for coUisionalities higher than the small value used here, many fewer modes in 9 and ^ 
are required, and execution time is then well below 1 second on a laptop. 

Further discussion of code benchmarking and results can be found in Ref. 

5. Multiple species 

In computing the bootstrap current, the intcr-specics coupling was formally negli- 
gible, but for other applications it is essential to solve kinetic equations for multiple 
species simultaneously, retaining cross-species interaction. The thermal speeds of differ- 
ent species may be quite different, so a grid in v which is appropriate for one species may 
not be capable of resolving another species. To maintain the high accuracy discussed 
in Section [21 the distribution function of each species fai should be represented on the 
zeros of the Pn{xa) using the normalized speed Xa for that particular species, as shown in 
figure [TT] Then the field term in the collision operator will require mapping between 
the Xa grid and grid and vice- versa. (The test particle term only involves /f,o , not /hi , 
so no remapping of unknown functions is required, and the test particle operator leads 
to a matrix that is diagonal in species.) Let us consider how to accurately perform these 
mappings. 

The last term in p2p . involving /^i and not the potentials, requires mapping the 

Maxwellian-like function /hi from the Xh grid to the Xa grid. This corresponds to mapping 

between grids 1 and 3 in fi gur e [TTl This mapping may be accurately done using the polint 

algorithm of Refs. 2lJ, [22, |23| , as discussed earlier for mapping from the polynomial grid 

to the uniform grid. This algorithm works very well both for evaluating values of fbi 
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• Grid 4: for H and G 
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Figure 11: (Color online) For calculations involving S species, 25 grids in the dimensional speed v are 
required: one for each distribution, plus one for the Rosenbluth potentials of each distribution. The 
polynomial grids for the distributions and the uniform grids for the potentials are each scaled to the 
thermal speed of the associated species. Hero, grids arc shown for a slow species (subscript s) and a fast 
species (subscript /). 



27 



close to Xb = (which is necessary when Vf, 3> Va, mapping a function known on grid 
1 to grid 3) and also for extrapolating to large values of Xb (which is necessary when 
Vb ^ Va, mapping a function known on grid 3 to grid 1.) Extrapolation is accurate 
because functional behavior cxp(— a;^) for large Xb is naturally incorporated into the 
polint algorithm. 

Now consider the Rosenbluth potential terms in p^ . When Vb ^ Va, the potentials 
must be evaluated close to Xb = 0. This corresponds to mapping a function known on 
grid 2 to grid 3 in figure 1111 Cubic interpolation is effective for this mapping, just as 
in the single-species case. However, when Vb ^ Va, we must extrapolate the potentials 
to evaluate them beyond the grid on which they are represented. This corresponds to 
mapping a function known on grid 4 to grid 1. As discussed in section [31 the H potential 
behaves as oc l/x^^^ for Xb — > oo. We find in practice it can be very inaccurate to use the 
simplistic extrapolation H = for Xb > xmux- (In particular, the code described in the 
following section does not converge well if this extrapolation method is used.) However, 
the more accurate extrapolation Hi{xb) ~ {xuayi/ XbY^^ Hi{xMiix) for Xb > XMax works 
well and leads to good convergence, as we will demonstrate in the next section. 

When Vb ^ Va, we must also extrapolate (PG/dxl to form the first left-hand-side 
term in (fTS]) . Notice from Section [3] that for large Xb, each Legendre mode £ behaves like 

+ (26) 

for constants Kg and Lf . It is tempting to make the further approximation (PGe/dx^ 
Lex^^~^, suggesting the extrapolation tPGi/dxl = {xMax/xbY^^ [d'^Gi/dxfj^ for 
Xb > xmux- However, since the Ki term is only algebraically small compared to the 
term, this extrapolation requires a much higher XMax (and therefore many more points 
in the uniform grid for the potentials) than a more careful extrapolation which accounts 
for both the Kg and Lg terms. To derive such an extrapolation, we write (|26p at the the 
last and penultimate points in the uniform grid for the potentials, which we denote by 
Xbi and Xbj respectively. Eliminating Ki and Lg by these two equations, we obtain the 
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improved extrapolation 





(27) 



Wc find that this result, which we use for calculations presented below, leads to conver- 
gence for much lower XMax (^4 — 5). 

Due to the factor fao in ([T^ , every term in ([T^ will be exponentially small for large 
Xa, even if Va- Therefore the polynomial Xa grid is able to accurately resolve the 

field particle operator when vi, ^ Va- However, if Vb ^ Wa, features in f(,i will be mapped 
to small Xa, leading to sharp structures that may not be well resolved by the Xa grid. For 
this reason, more grid points in Xa may be required when the thermal speeds are very 
different from each other. Fortunately, as shown in Figure O the new polynomials have 
zeros clustered near x ~ 0, so grid resolution is quite good in this region, much better 
than with a Laguerre or associated Laguerre grid. As we will show in the following 
section, a demonstration for species with very disparate thermal speeds shows that the 
number of grid points in x required for good convergence remains quite small. 

6. Application 2: Impurity flow 

As an example application involving multiple species, we consider the neoclassical 
flows of the main ions and a non-trace impurity species in a tokamak plasma. Due to 
the low mass of electrons, electrons do not significantly affect the motion of the ions, so 
we do not need to include electrons in the calculation. The perpendicular flow and the 
variation of the parallel flow on each flux surface (the so-called Pfirsch-Schliiter flow) can 
be written explicitly, so it is only the average parallel flow on a flux surface that must be 
determined from kinetic theory. Following convention, the specific form of the average 
we compute is {V ■ B). 

For this example calculation, we consider a plasma consisting of deuterium and 
Mo+''^°. We choose the latter since it is often a plentiful impurity species in metal- 
wall tokamaks, and its large mass (jtimo = 48mD = 96 a.m.u.) tests the perfor- 
mance of the numerical method when species have disparate thermal speeds. We con- 
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sider an impurity density tt-mo/'^d = 0.0012, corresponding to Z^s = 2, where ZcS = 
(tid + ^Mo"'Mo)/("-D + Zuonuo)- Since Zcs — 1 is not <C 1, the impurities are not in the 
trace limit, and so they will have a significant effect on the main-ion flow. 

Contributions to the parallel flows arise from radial gradients in the density and 
temperature of each species, as well as from the radial electric field. For this example, we 
will compute the contribution from the density gradients, i.e. we formally set dTiy/dip = 
0, dTuo/dip = 0, and d^/dijj = 0. (Contributions from these other drives could be 
superposed later due to the linearity of the kinetic equations.) We assume the density 
scale lengths are equal: n^) / (driD / dip) ~ umo / {driMo / dip) . We also take the temperatures 
of both species to be equal (Td = Tmo = T) . 

We normalize both flows with the same factor: 

{Vu-B) = t/,^^ (28) 
eTiD dip 

(Vmo-B) = C/mo — ^ (29) 
eriY) dip 

and so our goal is to compute the unknown dimensionless quantities Ud and Umo- To do 
so, we flrst solve the appropriate normalized drift-kinetic equations 

2 -x^- 1 dB 

= — T?^P D 

2 ° 52 de 



and 



^ dfuoi i'^~^'^)xMo I dB \ dfuoi f' Cmo , . 

""^"^ 86 2B [de) qR^b-Ve ^^^^ 

Q Mo ^~ ^ 2 2;-[L^ 



ZMom2 no 2 ""^^^ B2 dQ 



for the distribution functions /di and /moI, where 



and 



Cd = Cdd{/di, /do} + CddI/do, /di} (32) 

+ CdMo{/di, /moo} + CdMoI/dO: /moi}, 

Cmo = CmoMoI/moi, /moo} + C'moMo{./moo, /moi} (33) 

+ C'mod{/mo1, /do} + CmodI/moO, /di}- 
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Figure 12: (Color online) Convergence of Uu and Umo - the normalized parallel flows of and Mo"'"^'' 
- in the multiple-species drift-kinetic code with respect to numerical parameters defined as in Figure 
1101 Except for quantities being scanned, the numerical parameters are Nx = 8, Ny = 80, A'^^ = 40, 
NL = 4, N8 = 10, and x^ax = 5. Good convergence at very small values of Nx is again apparent. The 
coUisionality is = 0.01. 

This time, = vj^qR^/vY) in both equations pO^ -pi p . where vd = 4'\/27rnDe'* In A/(3^mDr 
FinaUy, we compute the moments 

C^D = -^(b dxjyxl I ^dU.hi'^ (34) 

Umo = ( —^(-8/ dxMoXMo [ d^^/MoiV (35) 

ttV^ \mMoJ nuo \ Jo J-i I 

We again use the MiUer geometry detailed in the appendix. A preconditioning matrix 
is employed, obtained by dropping not only the blocks that are off-diagonal in x, as in 
the single-species code, but also dropping all blocks that are off-diagonal in species. 

Figure [T^] shows the convergence of the two-species calculation. Calculations are 
performed for z/*d = v' J'c'l'^ ~ 0.01. For simplicity, in each run we use the same resolution 
for both species. As for the single-species case, convergence is achieved at a very small 
number of grid points in x, with 1% accuracy achieved by Nx = 5. Higher Ny is required, 
but this causes no decrease in code performance, as shown in the figure, for the reasons 
discussed in Section [31 

31 



7. Discussion 



A common feature of kinetic simulations is the need to accurately evaluate both 
integrals and derivatives of Maxwellian-like functions using a single grid or modal dis- 
cretizaton. In this work we have presented a set of methods which allow very accurate 
integration and differentiation for such functions even with very coarse grids. The central 
concept is to employ polynomial approximation using the weight exp(— x^) and interval 
[0,cxj), which is equivalent to an expansion in the new orthogonal polynomials JH-JS). 
The Stieltjes procedure (detailed in Refs. Q, Q Q Q or Section 4.6.2-3 of Ref. Q, 
and summarized in Section [3] of the present paper) is used to generate the grid and in- 
tegration weights. Then the poldif and polint routines of Refs. 
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are used to 



generate the differentiation and interpolation matrices. We emphasize that these algo- 
rithms from previous references are applicable to any weight function and any integration 
interval, and were not developed with velocity-space applications in mind. The recog- 
nition that these algorithms may be applied to kinetic calculations, using the interval 
and weight appropriate for the v or vj_ coordinate of Maxwellian-like functions, is novel. 
Figure is a central result of our work, showing the new method performs far better than 
other discretization schemes at both integrating and differentiating functions relevant to 
kinetic theory. The abscissae, integration weights, and differentiation matrices are not 
significantly more difficult to compute than those for classical orthogonal polynomials, 
and public-domain source code is available for generating the grid, integration weights, 
differentiation matrices, and interpolation matrices. 

We have also shown how the discretization scheme may be applied when the Fokker- 
Planck collision operator is included in kinetic calculations. The Roscnbluth potentials 
in the collision operator cannot be discretized in the same manner as the distribution 
hmction because they do not behave as oc exp(— x^) for large x. Instead, they may be 
discretized on a uniform or Chebyshev grid with very high resolution. By the construction 
(|18p , this high resolution grid is eliminated through matrix multiplication when the final 
matrix for the kinetic equation is constructed, and so the grid resolution for the potentials 
does not affect the size of the final matrix. Consequently, overall code performance and 
memory requirements are essentially independent of the resolution of the grid for the 

potentials, and the final matrix size is determined instead by the much coarser resolution 
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of the grid for the distribution function. 

The methods have been demonstrated in several neoclassical calculations for a toka- 
mak plasma, which involve the solution of one or more drift-kinetic equations including 
collisions. Figure [TUl illustrates such a calculation for the bootstrap current, showing that 
as few as four grid points in x may be sufficient for convergence of two significant digits 
(which is more than enough accuracy for comparison to experiment.) When species with 
disparate masses are included simultaneously, as in the impurity flow computation of 
figure VTA careful extrapolation of the Rosenbluth potentials must be performed, but as 
long as this is done, convergence may be achieved for a similar number of grid points in 

X. 



Acknowledgements 

The authors are grateful to Michael Barnes for assistance regarding the GS2 grid. This 
work was supported by the Fusion Energy Postdoctoral Research Program administered 
by the Oak Ridge Institute for Science and Education. 



Appendix A. Miller geometry 

Here we present several relations for Miller equilibriajincluding some useful relations 
that are not written explicitly in the original reference |3l| . Our focus is to compute the 
quantities needed for the neoclassical calculations described in the main text, B - VO and 
the field magnitude S, both as functions of a poloidal angle 9. 

The Miller formulation begins with the assumption that the flux surface of interest 
has the shape 

R{e) = Ro + rcosie + xsme), (A.l) 
Z{e) = nrsine, (A.2) 

where Rq is the major radius at the geometric center of the flux surface, r is the minor 
radius, k is the elongation, x = arcsinJ, and S is the triangularity. The aspect ratio 
A may be defined by ^ = Ro/r. All the quantities Rq, r, k, 6, and A are considered 
to be functions of the fiux surface label ip. This 9 angle is generally not the angle 
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arctan(Z/(i? — Rq)), nor is it the Boozer or Hamada angle or other straight-field-hne 
coordinate. 

Next, consider the standard relation for any {q^,q^,q'^) curvilinear coordinates 
where 

^-fr'^xf-W v,^xv,»)- (A,4) 

is the Jacobian, and here we choose (g^, g^, g'^) — (i-'X^f^) with ( the toroidal angle. 
Recalling that any axisymmetric field may be written B = VC x Vip + ^VC, the quantity 
B ■ W9 we are seeking is equal to 1/J. To evaluate the derivatives in (|A.3|) - (|A.4p . we 
write the position vector as 

X = exRcosC, + eYRsm( + ezZ, (A. 5) 

where (ex, ey, ez) are the Cartesian unit vectors. Differentiating (|A.5[) . we obtain, after 
some algebra, 

1 Y 

" B~V9 = ^^-^^ 

where 

Y = Kri? |sin(6' + a; sin 6*) [f + a: cos 6*] (f + s^) sin 6* (A. 7) 

" OR 

— h cos(6' + X sin 9) — ss sm{9 + x sin 0) sin ( 

or 

Sk — {r / H)dK/ dr and ss = {1 — 6^)^^/^r dS /dr. Then from the magnitude of (jA.3[) . 

IVV'I = ^\/{sin(6i + a;sin6i) [1 + x cos 61]}^ + (k cos 61)2. (A.8) 
Defining Bq — I/Rq, the field magnitude may be written 



B = JB^g + (BoRo/R)^ (A.9) 



where Be = |V?/'|/i? = (jA.8|) /i? is the poloidal field. The Bg calculated in this manner 

from (|X6l) -(|X8 l) is precisely Eq. (37) of Ref. [sil- 

Equations (|A.6p - (|A.9p with (|A.ip are the desired expressions for the two quantities 

we need, B and B ■ V9, albeit in terms of the quantity dijj/dr in (jA.6p . It is usually 
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preferable to specify the safety factor q in place of dip /dr. To do so, we observe 



27r L de 2n n B - Ve ^ ^ 



Then using (jA.6p . 

dip BqRo f^"" ^^Y 



, dO— (A.ll) 

where the integral may be evaluated numerically. 

In summary, the magnetic geometry is fully specified by the following seven dimcn- 
sionless parameters: A (aspect ratio), k (elongation), 5 (triangularity), q (safety factor), 
Sk (elongation shear), ss (triangularity shear), and dRo/dr (related to the Shafranov 
shift). Since only the normalized quantities B = B/Bq and qR^b ■ VO appear in the 
code, there is no need to specify the two dimensional parameters r (minor radius) and 
-Bo (magnetic field at the geometric center of the fiux surface.) The additional dimcn- 
sionless parameters s (global magnetic shear) and a (normalized pressure gradient) are 
discussed in Ref. jsij but are not needed for the calculations here. 

The calculations for figures [TOl and [T2] were performed using the following parameters. 



identical to Ref. 



311: A = 3.17, K = 1.66,(5 = 0.416, s„ = 0.70, = 1.37, dRo/dr 



-0.354, and q = 3.03. 



References 

[1] F. Valentini, P. Veltri, A. Mangeney, J. Comp. Phys. 210 (2005) 730-751. 

[2] M. Barnes, W. Dorland, T. Tatsuno, Phys. Plasmas 17 (2010) 032106. 

[3] R. Numata, G. G. Howes, T. Tatsuno, M. Barnes, W. Dorland, J. Comp. Phys. 229 (2010) 9347- 
9372. 

[4] A. Pataki, L. Greengard, J. Comp. Phys. 230 (2011) 7840-7852. 

[5] E. A. Belli, J. Candy, Plasma Phys. Controlled Fusion 54 (2012) 015015. 

[6] C. W. Clenshaw, A. R. Curtis, Numerische Mathematik 2 (1960) 197-205. 

[7] L. N. Trofethen, Spectral methods in Matlab, SIAM, Philadelphia, 2000. 

[8] M. N. Rosenbluth, W. M. MacDonald, D. L. Judd, Phys. Rev. 107 (1957) 1-6. 

[9] L. Pareschi, G. Russo, G. Toscani, J. Comp. Phys. 165 (2000) 216-236. 

[10] Z. Xiong, R. H. Cohen, T. D. Rognlien, X. Q. Xu, J. Comp. Phys. 227 (2008) 7192-7205. 

[11] J. A. Spencer, J.-Y. Ji, E. D. Held, J. Comp. Phys. 231 (2012) 6192-6206. 

[12] S. K. Wong, V. S. Chan, Plasma Phys. Controlled Fusion 53 (2011) 095005. 

[13] B. C. Lyons, S. C. Jardin, J. J. Ramos, Phys. Plasmas 19 (2012) 082515. 

35 



[14] M. Barnes, I. G. Abel, W. Dorland, D. R. Ernst, G. W. Hammett, P. Ricci, B. N. Rogers, A. A. 

Schckochihin, T. Tatsuno, Phys. Plasmas 16 (2009) 072107. 
[15] W. Gautschi, Math. Comp. 22 (1968) 251-270. 
[16] W. Gautschi, SIAM J. Sci. Stat. Comp. 2 (1982) 289-317. 
[17] W. Gautschi, ACM Trans. Math. Software 20 (1994) 21-62. 

[18] W. Gautschi, ORTHPOL, Accessed October 6, 2012. http://www.netlib.org/toms/726. 

[19] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Webnote No. 3, Rev. 1: Implemen- 
tation of Stiel, Accessed October 6, 2012. http://www.nr.com/wcbnotes73. 

[20] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical recipes, 3rd. cd., Cam- 
bridge University Press, 2007. 

[21] J. A. C. Weideman, S. C. Reddy, ACM Trans. Math. Software 26 (2000) 465-519. 

[22] J. A. C. Weideman, S. C. Reddy, A Matlab Differentiation Matrix Suite, Accessed October 6, 2012. 
http:/ /dip.sun.ac.za/ wcidoman/rcscarch/diffcr.html. 

[23] J. A. C. Weideman, S. C. Reddy, DMSUITE: Matlab Differentiation Matrix Suite, Accessed October 
6, 2012. http://www.mathworks.com/matlabcentral/filecxchange/29. 

[24] P. Helander, D. J. Sigmar, CoUisional Transport in Magnetized Plasmas, Cambridge University 
Press, Cambridge, 2002. 

[25] B. Li, D. R. Ernst, Phys. Rev. Lett. 106 (2011) 195002. 

[26] P. J. Catto, K. T. Tsang, Phys. Fluids 20 (1977) 396-401. 

[27] J. G. Cordey, Nucl. Fusion 16 (1976) 499-507. 

[28] J. B. Parker, P. J. Catto, Plasma Phys. Controlled Fusion 54 (2012) 085011. 

[29] F. L. Hinton, R. D. Hazeltine, Rev. Mod. Phys. 48 (1976) 239-308. 

[30] O. Sauter, C. Angioni, Y. R. Lin-Liu, Phys. Plasmas 6 (1999) 2834-2839. 

[31] R. L. Miller, M. S. Chu, J. M. Greene, Y. R. Lin-Liu, R. E. Waltz, Phys. Plasmas 5 (1998) 973-978. 
[32] G. L. G. Sleijpen, D. R. Fokkema, Electr. Trans. Num. Anal. 1 (1993) 11-32. 

[33] S. Balay, J. Brown, K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Kneplcy, L. C. Mclnnes, B. F. 

Smith, H. Zhang, PETSc Web page. Accessed October 6, 2012. http://www.mcs.anl.gov/petsc. 
[34] S. Balay, J. Brown, , K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. 

Mclnnes, B. F. Smith, H. Zhang, PETSc Users Manual, Technical Report ANL-95/11 - Revision 

3.3, Argonne National Laboratory, 2012. 
[35] F. L. Hinton, M. N. Rosenbluth, Phys. Fluids 16 (1973) 836-854. 
[36] M. Landreman, D. R. Ernst, Plasma Phys. Controlled Fusion 54 (2012) 115006. 



36 



